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= cross-sectional area for main steam floA^r 
and secondary steam flow at the throttle 
valves (in^O 

= constant which relates the flow aid pressure 
drop in a tLirbine with many stages 
= coefficient used in determining steam flow 
through main and secondary throttle valves 
('l\''sec) 

= heat transfer coefficient of feed water 
heaters (ws/lbl 

= heat transfer coefficient of reteater 
(ws/lb°F) 




^f» ^fw’^fw’^0 “ 


J, Kj 
ICggp, K^p 


enthalpy of main steam at nozzle chest, 
reheater, throttle valve, moist-ure separator 
and isentropic end points from pressure 
and respectively (wsAb) 
enthalpy of saturated water in the reheater 
feed water leaving heater 1 , heater 2 and 
water in the condenser r espectiv^y(wsA'K)) 
constant to relate mechaiical energy and 
vice versa (ws/ftlb,ftlb/ws) 
coefficient ^ich indicates the amount of 
bled steam from the H.P. and L.P, turbines, 
respectively. 
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constant in the Callender’ s empirical 

equaticn relating pressure, density and 
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enthalpy of superheated steam (Ih^ft / 
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reheater (ib/in.^) 

heat transferred in heater no. 1 and 2, 
respectively (vs) 
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time ccnstants associated with heater no. 1, 
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torque of high press-ure and low pressure 

turbines respectively (lb ft) 
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total effective volume of the nozzle 
chest and reheater, respectively (ft^) 
specific volume of saturated water and 
steam, respectively, around a pressure 
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, W 2 , = steajtn flow rate at tlirottle .valve, 

nozzle chest, L«P« turbine, before 
the reheater and after the reheater 
respectively (Ib^^/sec) 

WfiKp? = bled steam frcm high pressure and 
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= efficiency correction factors to take 
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J^R, /"a 
J1 


density of steam in the nozzle chest, 
in the reheater and at the turbine 
outlet (Ib/ft^) 

speed of turbine equal to 120 (rad/seo). 



SYNOPSIS 


A method is developed for low order modelling 
of large maltivariable systems. The method is applicable 
to both linear as well as non-linear s^'-stems and does 
not require the transfer f inction or vector differential 
equation of the system to be specified. Onlv measured 
input-output data cf the system at discrete intervals 
of time is required to determine its linear^ low order, 
discrete- time model. The method is based cn the use of 
matrix pseudo-inverse to estimate the parameters of the 
model which minimize the sum of the squares of the errors 
between the responses of the original system and the 
low order model at the sampling instants. A recursive 
algor ithim is used which reduces the computational and 
data storage requirements and enables successively 
improved model of the system to be determhied as new 
data from the system is made available. Thus the method 
can be used for on-line application s. 

■ki iterative technique is used to convert the 
discrete-time model into a corresponding continuous- 
time model* 

The method is used to obtain low order model of 
a nuclear reactor turbine system- 



CHjE>TSl I 


1 «1 In u-^^o d- uction 

Modern systems are often large, resiiLting in 
dynamic model which consists of a large number of non- 
linear differential equations* It is common practice to 
linearize the model around an operating point and assume a 
time-invariant behavior# The resulting linear model describ 
the behavior of the system adequately in the neighborhood 
of the operating point and in terms of computational effort 
the system simplifies considerably* In many situations, 
such as a thermal power plant, a nuclear reactor system, 
a process plant or aircraft systems etc*, the dynamic 
model is still of such a large order as to be inconvenient 
or impractical for its simulation or control system design* 
Thus, it sometimes becomes important to look for a dynami- 
cally equivalent model which is of reduced order. The 
reduced order model should describe the behavior of the 
system almost as well as mere complicated higher order 
model and retain the most important parameters of the 
original system* 

In the last decade or so considerable attention 
has been devoted to the problem of deriving reduced- order 
models for complex systems. The problem is clearly of 
great interest to those working in the applications area 
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1 *2 Uses of reduced order models 

ifeduced- order models are obtained for the 
folio-wing reasons : 


(1 ) To simplify understanding of the system, whether the 
problem is analysis, synthesis, or identification and to 
achieve easy simulation of the processes# 


(2) To reduce computational requirements and to make 
better use of the supporting facilities as well as the 
ccmputer itself, in addition to widening the scope of the 
problems which can be handled on a given size of machine. 
Of particular significance to on-line computer control 
applications is the reduction of the computational effcxrt 
of optimal and adaptive controllers by deriving suboptimal 
strategies based on reduced order models\__4-J, 1^53 • 


(3) To obtain a lower dimensional ’’control law” for 
simplifying the structure of the feedback controller l6l. 


0 +) To make best use of scanty experimental data, by 
estimating a few parameter with confidence, rather than 
estimating more parameters with less confidence. The 
low order model is then a more reliable predictor Ln- 

(5) To generalize results established on a particular 
system to comparable systems, using low order models 
which basically differ in time scale only. This is 
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paJ'ticularly usefiil in tolerancing performance for system 

8 -}. 

1 »3 Liter at ure survey 

Although most of the work in the area of low 
order modelling of complex systems has been done in the 
last decade or so, the subject has aroused widespread 
interest. In the beginning different investigators 
presented radically new approaches to the same problem. 
Later on several variations and combinations of these 
approaches were developed. These different approaches 
along with their modifications are presented below grouped 
under different categories* 

1 .3.1 Mpdelling from loop data 

The initial Impetus came from Axelby ^ 9*^ and 
Biernson 0j[ who developed methods for rapidly estimating 
system poles from the open loop Bode plot. This approach 
was formalised by Kan Chen who suggested that 

dipoles might be estimated from system behavior in the 
high gain region, fax-off poles from the low gain region, 
and dominant poles from the cross ovcx region. Kan Chen 
method has been' exploited by CanflGld{[i2']for multiloop 
systems, and by Towill ^133 in terms of a rational design 
procedure* The basic problem with these methods is that 


check-out 



Ij. 

if there are many break points in the high and low gain 
regions, the dominant poles, and hence the transient 
response of the low order model, is considerably in 
errcr* 

1 «'^«2 Modelling from the g y stem tran s fer function 

Since the system transfer function contains all 
the data needed to model the system, there has always been 
effort to use the coefficients directly in a low order 
model* Deriving the system transfejr function coefficients 
is time consuming if a computer sub-routine is not available, 
so that for high order systems it is not a preferred 
approach for hand computations* Oldenburger neglected 

seme of the coefficients on an intuitive partitioning 
basis and applied it to develop a reduced order model of 
a turbine governer system. Sil 3 ak\_ 15 '^ , in parameter 
plane design methods, simply neglected the higher oid.er 
coefficients and varied the lower order coefficients by 
adding a compensator to yield the original system response. 
The method due to Chen and Shieh is based on continued 

fraction expansion (GFE) of the transfer function of the 
original S 3 ^sfcem* The numerator and denorainatcr pol37nomials 
in the transfer function are arranged in the ascending 
order, in powers of s, and then expanded into the Cauer- 
type continued fraction. From a consideration of the 
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final valuG theorem, it follows that quotients in the 

expansion are in the order of decreasing significance* 

Thus truncating the CFE at a suitable stage gives the 

reduced transfer function of the desired accuracy. The 

GFB used is equivalent to Taylor's series expansicn about 

s = 0 and can be shown to be a special form of Fade' 

approximatiai for a scalor rational function 1.171. The 

method ensures that the model gives the correct steady 

i 

state response for a class of inputs, namely^ t , 
i = 0, 1, 2,*.*p where 'i* is a positive integer depen- 
ding on the or'der of the model. However the initial 
transient response may not be good. Shamash^l , 
expanded the numerator and denominator into a Jordan-type 
continued fraction, about s = 0. This method requires 
much less computation than the original method. Ghuang*^9^-, 
1203 modified Chen and Shieh' s method. He expanded the 
system transfer function into a Gauer-type continued 
fraction about s = 0 and s = 00 . By alternately expanding 
around s = 0 and 3 = OO , ‘this method ensures equally 
good steady state and initial transient responses. 

However, biased responses m^ also be obtained. By 
'fitting’ more coefficients in the power series expan- 
sion about s = 0 or s = OO , a better steadir state or 
initial transient response respectively, may be obtained. 
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Ghen t2G extended his method to multivariahle systems 
by expressing the transfer function matrix as a fraction 
with a matrix polynomial ic the numerator and a scalor 
polynomial in the denominator in the ascending order and 
then expanding into matrix continued fraction as in the 
case of scalor systems* Truncating the continued fraction 
at a suitable s tage gives the red-iced order multivariable 
model. Chen applied this method to reduce the model of a 
gas-turbine system. 

Vittal fiao and Lamba minimised the integral 

of the weighted error function between a given transf^ 
fmction and the reduced transfer function in the frequency 
domain. This method is of especial use in the non-linear 
systems where the frequency response of the linear part 
of the simplified system should approximate the frequency 
response of the linear part of the original system. 

Shamash l23] extended the continued fraction 
methods for linear continuous- time systems to linear 
discrete- time systems. 


1 «3.3 Eigenvalue apuroach 

Evans is considered the originator of this 

approach as his root locus method brought the determination 
of eigenvalues into the realms of practicability for 
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realistic order systems* The first computerised, state- 
space approach to the problem is credited to Davison 
and essentially consists of neglecting eigenvalues of the 
original system farthest from the origin, retaining only 
the dominant time caistants* The retention of dominant 
eigenvalues makes the response of the reduced model 
approximate that of bhe original, since the eigenvalues 
neglected make a very insignificant contribution to the 
total response except at the beginning. Relationships 
frcxn the time-solution of the vector differential equation 
of the original model were used to develop a reduced model 
which maintained both the correct proportion of the eigen- 
vectors and the desired eigenvalues* The resulting low 
order model did not preserve steady state values, although 
the initial and middle parts of the responses were good 
representations. Marshall hy a similar reduction 

procedure and by neglecting the dynamics of fast modes 
of the originsflL high order model, derived a lower order 
model which additionally achieved steady state compara- 
bility also and offered zero steady state error for cons- 
tant inputs. Following the publication of Daviscn's work, 
there was much correspondence between Davison and 
Chidambara , but according to Marshall “^27*} , no new 
method was proposed and the correspondence consisted 
essentially of variations upon the known theme. 
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Chid ambar a 128] accounted for the steady state error in 
Davison’s method by pointing out the terms which cannot 
be neglected as steady sta.te is approached. However, 
Davison 283 pointed out that Chid ambar a’ s method while 
giving correct steady state values, did not give good 
transient response agreement. This was because Chidambara 
was essentially finding a reduced forcing function and in 
doing so it did not excite modes of t’ e approximate model 
in the same proportions as found in the original model* 
Chidambara I 293 •, 1301 ccf.;ntered Davison ' s arguments by 
giving exanples* However, Davison^29'^ , ^ 3 0][ pointed out 
that the examples given by Chidambara were misleading 
because only one mode was presented in the simplified 
model. He also mentioned that Chidambara* s method was 
more difficult to use in practice. Following these 
publications, Davison ^31 j modified his method by defining 
a diagonal matrix and incorporating it in his original 
method in sudi a manner that the system gave satisfactory 
dynamic response and correct steady-state response for 
step inputs. Inspite of a number of these articles, it 
was not clear which model should be used since the examples 
given showed that either method is correct or incorrect 
according to various cases considered, FossardY.32^ noted 
this and pointed out that the problem lied in the existence 
of null eigenvalues. He modified Davisons method to tahe 
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care of null eigenv^ues and produced satisfactory steady 
state and dynaniic responses for step input. He extended 
his method to general input case also by observing that 
if the slow modes of the system are correctly excited, 
those of the input will be incorrectly transmitted unless 
they are integrated into the system and the reduction 
process made on the augmented system. Newell and Fishei^S] 
applied Davison’s method to reduce the order of an evapora- 
tor model. Seborg and Fisher extended the continuous- 
time methods of Davison and Marshall to discrete timo 

systems. Jn interesting variation to Davison's method is 
that proposed by Kuppurajulu and ELangovan ^53 'wherein 
the high order system is replaced by three models, 
successively representing the initial, intermediate and 
final stages of transient response by identifying those 
eigenvalues which are most dominant in the initial, inter- 
mediate and final portions of the system response. 

Mitra' s 31363 method consisted of orthogonal projection 
onto a subspace of the state space determined by a defined 
matrix. M.lscn 3,373 determined a low order model by 
minimising a functional of error between the output vector 
of the system and the output vector of the model. The 
algorithm is based on the solution of a Lyapunov 
matrix equation and requires optimization of individual 
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parameters occurring in the particular canonical .fprm 
chosen* However this method is cumbersome and quite time 
consuming. He extended it to multivariable systems also 
W and used it on a boiler model. 

1 *3«H- Modelling from response data 

This is the general identification problem if we 
probe to find low order models which acceptably describe 
the ‘black box* under test, through its input-output data. 

In the frequency domain, the historical path commenced 
with weighted least squares method of Levy \39j • He 
fitted the frequency dependent polynomial into the generated 
data by minimising an error function whic±i was a measure 
of error in the fit. The method has been widely used and 
is capable of fitting some quite unusual frequency 
responses , but several difficulities have been noted 

with the method. Sanathanai and Koerner^l"] found that 
the procedure described does not give a good fit if the 
frequency data spans several decades and have proposed 
instead an iterative procedure. The authors cited an 
example in which the it eratively- computed parameters 
differ considerably from those obtained by Levy’s method. 
Sumner has made the same criticism of the Levy method, 

but proposed a technique in which the error between the 
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data and the model is minimised using a non-linear optimi- 
sation method. A further mcd.ification made to the method 
by Payne is to incorporate certain constraints based 

on additional toowledge of the system, such as steady 
state gain* He noted that in the unconstrained system 
there is a tendaicy for poles to occur on the right half 
of the complex plane and so give unstable responses in 
system Inoi^ to be stable. He reported, in fact, that in 
fitting several hundred responses about 30fo resulted in 
unstable systems when no constraints were used, but that 
the use of constraints reduced this to only 1 per cent* 
Towill and Payne ^44-1 studied many systems and came to the 
conclusion that the least squares error criterion did not 
necessarily result in a satisfactory model in many systems* 
However he found this criterion to be most straight- 
forward in terms of application and least computation time 
consuming in most of the cas s* Considerable time may 
therefore be spent in system study and simulation until 
a satisfactory compromise is reached, 

A number of methods minimising some function of 
the error between the system and the low order model at 
sampled response times have been developed, Aiderscn^V^’} 
determined a reduced order model, the response of which 
approaches that of the original in such a manner that the 
mean-square error between the two responses, over a given 
finite internal, is minimised* This is '’ccsomnli shed hv 
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using the orthogonal projection theorem in the theory of 
linear ■vector spaces to nrLniuiize the sum of the squares of 
errors betitfeen the responses of the actual system and a 
discretized model at the sampling instants* The vector 
dif ferentiaO. equation is then obtained from the vector 
difference equation of this discrete model. Sinha|%6*3 
derived a reduced order model by expressing the m.odel of 
desired order in terms of a pulse transfer function, with 
its coefficients as unlnowns. To find out these unknowns, 
he expressed the transfer function in an equivalent 
difference equation form. This equation when 'written for 
different time instances gave a set of ordinary algebraic 
equations with the coefficients of the reduced transfer 
function as unknowns* The equations were solved for 
these unkno'wns by calculating Left Minimum Pseudo-Inverse 
b8] which was shown to give the unknown coefficients 
in such a manner that the sum of squares of the error 
between the original and the reduced model responses was 
minimized at the sampling instances* Ai iterative algo- 
rithm was used and thus the method required much less 
storage of data and much less confutation than Anderson’s 
method. Bereznal and Sinha 1.^91 Proposed another method 
based on the minimization of deviation between the 
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responses of the original system and the reduced order 
model in the direction that is everywhere perpendicular 
to the response of the system. He considered various 
error minimization criteria. However, he noted that 
analytical solution to the problem of error minimization 
is available only for the least squares error criterion. 

He suggested the pattern search techniques of Hooke and 
Jeeves ^50^ solve problems formulated in other errcr 
minimization criteria. Meier and luenberger 
minimised the mean square difference between system output 
and model output when a random process is used to drive 
both systems. Hli£ 523 obtained a low order model from 
impulse response estimates by computing the time moments 
of the impulse response. The moments were then, using 
series expansion theory, equated to the unknown transfer 
function, and the coefficients determined by elimination. 

As the model order increases, considerable difficulty is 
experienced with convergence of the impulse moments for 
realistic test data. Nevertheless, it is a powerful 
technique for use in heavily damped systems* rlggs and 
BigarJ2531 developed a general framework for optimum 
model reduction based oa impulse responses of the original 
and reduced single input - single output systems and gave 
example of a heat diffusion system. Galiana suggested 

a method of approximating linear time invariant 



miilti input - multi output systems by a lower order system 
through a weighted least squares minimization of a seal or 
integral functioi of the impulse response matrix over all 
positive time. A set of necessary conditions was derived 
in the form of a number of algebraic matrix equations, the 
solution of which yielded the unl-snovjn system parameters. 
Although he gave several examples, he did not apply it 
to any real system* -Applevich' s '^55'^ raethed on multi- 
variable linear time- in variant discrete- time system was 
to optimally approximate the given system with respect to 
the sum of the squared output errors of the impulse response. 
He also established existence and uniqueness of the 
solution and in this connection developed certain projec- 
tion theorems. Although he did not apply his method to a 
practical system, he gave an exanple. 

1 .3.5 Other approaches for low order modelling 

Shi eh and Wei Proposed a method which was 

a combinaticn of dominant eigenvalue and continued frac- 
tion approaches. In their method which was applicable to 
multi-input multi- output systems, they expressed the 
transfer function matrix as a product of two polynomials. 

The dominant eigenvalues were usc^ to formulate tie common 
deiominatcr polynomial of a reduced order model aad the 
numerator dynamics of the reduced order model were 



obtained by matrix- continued fraction approach* The 
method has the advantage of eigenvalue approach in the 
sense that reduced model is always stable* It also has 
the advantage that the computationally unattractive 
procedure of complicated linear transformation, matrix 
diagcnalizaticn , and steady state value matching, etc. 
involved in eigenvalue ^proach is replaced by a simple 
technique of continued fracticn. Arumugam and Rama Moorty 
■5'^used Schwarz caicnical form in simplifying dynamics of 
systems described by linear differential equations* The 
method does not require the computation of eigenvalues aid 
eigenvectors* But this method cannot be used for 
multivariable systems* Ravi chandr an used singular 

perturbation mothod to reduce the model of a nuclear 
reactor system. 

1 *f|- Present work 

In the present work a method has been suggested 
to determine low order model of large multivariable system. 
This method is developed following the guidelines of 
Sinha’sf ^ 6 } work. Sinha developed a low order modelling 
technique from input-output data of a linear or non-linear 
system, measured at discrete intervals of time. Low order 
model is determined by using matrix pseudo-inverse 
technique which minimizes the squares of error between 
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the mo’ g 1 and system responses at the sampling instants* 
However, this method was exclusively for single-input 
single- output systems* Moreover, this method was not 
applied to any real system. Sinha selected a fictitious 
linear system with its parameters chosQi in such a 
manner as co simplify the analysis. 

The present method is applicable to multivariable 
systems and cm be used for on-line applications* The 
method is applicable to both linear and non-linear systems 
and does not require system transfer function or vector 
differential equation to be specified. Only the input- 
output data of the system measured at discrete intervals 
of time is required to determine the low order discrete 
time model. As in Sinha' s method, matrix pseudo- inverse 
is used to minimize the squares of error between the 
responses of the low order model and the original system 
at sampling instants* 

An iterative algorithm is used which reduces the 
computation and storage of data. The discrete time model 
is convert©! into corresponding continuous time model by 
application of Smith's L611 method . 

The method is used to obtain low order model of 
a nuclear reactor turbine system. 
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Chapter U describes the mathematical develop- 
ment of the present method and the algorithm for obtaining 
the lower order model. In Chapter III, the model of the 
nuclear reactor turbine (which has been used here as an 
example) has been discuss®! and application of the present 
method has been developed. In Chapter IV, the computa- 
tioial results ccnsidering various orders of reductions 
have been givoa. Concluding remarks, limitations of the 
application of the method aid scope for further work 
has been discussed in Chapter V. 



ciiiiPTEr: II 


P RQBLM FORlIJl^ATIon MD HETHO D OF SOL U TION 

We aim at the determination of a low order 
model of a multivariable system from its input- output 
data obtained at discrete intervals of time* Such a case 
might arise in modelling of a large system by a low order 
model, v/hich has considerably smaller order dynamics, by 
using directly the measured input -output data* 

Once the algorithm to find lower order model 
from actual input-output data has been incorporated in 
the computer, it will give successively improved lower 
order model as density of the input-output from the plant 
increases* This can be used as an on-line process* 

2*1 Problem formula tion 

The low order model of a system may be represented 
by a block diagram shown in Fig* 2.1. For the ’p' inputs 
(4), T 2 C^),**.rp there may be ’S’ outputs 

C-] (-4), O 2 OS), ••• Cg (A) of the model. The relationship 
between the input and output variables is described by 
transfer function matrix H (A), which in effect describes 
the low order model. 

If the system is modelled as low order linear 
model, each input can be treated independently of each 
other. Conplete output of the system can then be obtained 





Figure 2*1 Ilulti-Input Multi-Output System 




by superposition, i.e. , outputs corre spending to' each 
input alaie are added togj^ther. Thus, the i-th output 
(^) is given by, 

P 

C^) ; i = 1 (2.1) 

0=1 ~ 

where r^ (<%) is the j-th input and H^-^^is the transfer 
function between i-th output and o-th input with all 
other inputs reduced to zero. For different values of i 
and 0, this describes the transfer function matrix with 
each H^i^J)as an element transfer function. Now, the mode- 
lling process reduces to finding coefficients of each 
of the transfer function such that the output variables 
of the low order model approximate those of the original 
system. 

2.1.1 Low order modelling for discrete system data 

Since in practice it is easier to make measure- 
ments at discrete intervals of time than to measure them 
continuously, the system is m.odelled as low order discrete 
model. The model will be derived so that its output 
approximate those of the original s"‘>''stem at the sampling 
instants. Such a model can be described by a digital 
element wherein the output is in the form of pulse 


sequences. 
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In a digital element, the difference e':iuatiai 
relating the input and output pulse sequences is, 

=liLt<^2)TV ... 

\l,n =KlC 

®-ki,o * ^a,i ^ i) ^,2 

"ia,m (2-25 


+• • • + 


where , 


th 


r^ is the 1 input variable* 


th 


Ct^ is the change in the k oiitput variable 


■l-T-l 

due to 1 input alone' 


^kl,1»**'’ \l,n ‘^kl,o’'**'=^kl,m 
various coefficients* 


and ar 


etc* are the 


T is the sampling interval 
i is an integer* 

This linear relation can be interpreted as a 
formula through which the present output number cj^'iT) 
can be computtai by taking weighted sums of a fixed group 
of input and output numbers* The process of digital 
element is more clearly depicted in the Fig* 2*2* 

The output c^2. t = iT depends on 

input q. -applied at t = iT as well as on the previous outputs 
at times (i~1)T, (i-2)T,**. (i-.n)T and on the previous 


# * « 
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inputs applif-d at tines (i-1)T, (i-.2)T,. , (i-m)!. The 
relative effects of these variables in the previous time 
values on the present output variable c^(iT) is given 
by the weighting coefficients • • *b^^^ and a^ q... 

^kl,m* 

In a linear multi-input multi- output system, 

superposition theorem can be applied; so that the outputs 

corresponding to each input alone are added together to 

give the total output* In other words, 

P 

0^ (IT) = Xl V (2,3) 

1=1 

where , 

(IT) is the lo-th total output at t = iT when 
’p’ inputs are applied* 

C]^ (iT) is the k-th output at t = iT correspen- 
«iing to 1-th input alone. 

So, in the multivariable, linear discrete system, the 

appropriate difference equation which relates the input 

and output pulse sequences is obtained by taking summaticn 

on both sides of eqn. (2.2) from 1 to p, and is, 

P p 

= H 0^ Cil] = 

1=1 1=1 ^ 

1=1 


+...+ bki^^ c^£(i-n)fi); k =:1,2,...,s 


(2.4) 
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2»1>2 Zp-transform analysis of the system 

3h order to describe the dynamics of the discrete 
process represented by the linear difference equations 
given in the previous section, Z- transform analysis is 
tised. 


Multiplying both sides of eqn. (2.2) by 
where s is the Laplace-transf orm variable, and taking the 
siimmatiai from i = 0 t o qo , we have , 

ifs ^ r., -iTs 


i=0 


CO 

i=0 


- 

i=0 

t~, -n -ifs 

T \l,n ^klO^-’") H 


• • ♦ — 




. r Cille"^^^ + 


kl,o 1’ 


i=0 

oo 


^ ®kl,1 ^it ® ^1 2 


1=0 


-ITs 


i=0 


i=0 


pO 

*—* Y . ^KL,b <2*« 

i=0 

The lower ir;dex in the summation of (2.5^) is zero, 
which indicates the fact that the values of the sequences 
for the negative time are zero. 
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Putting in eqn* (^.* 5)5 Z = e Tuom the defi- 
nition of Z- transform and maiding use of the 'Shifting 
theorem' which states, 




and Z"^ 


00 


00 


51 - >^^,1 H. Z 


-i-1 


i=0 

— • • • — h 


oO 



Z— 

i=0 

CsO 


firlz 


i=0 


00 

i=0 


.i-2 


CO 


a 




■f 




■i-1 


hl.o / 1 

^ tier- 44b 

i=o 

CJO 

+,»»+a, T* rXii^Z-^ 


-m 


i=0 


\;l,m I — ^”1' 
i=0 


( 2 . 6 ) 


Factoring out the common summation for each of the terms 
of eqn* (2.6), these result'’, 

- -i -1 -2 

i=0 


00 

n 


^^1,0 

^^,1 z'”' -’-••• + a^,„.Z““5 

i=0 





Again from the 

definition of Z- transform. 


c_ (Z) 
Id. 





i=0 




ttO 


and 

(Z) 

= i: 




0 

11 



(2o7) 
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Using '"Uese definitions in eqn® (2o7)j have, 


c^^-,(Z) (1-h. 


hl,1 


Z“'l - b. 


klj2 


Z 


V,n 2‘U = 


■1 — m 


m 


or, 


(2) a^-i 


o 


+ a' 


,-1 




,-m 


ia,1 -K..+ a^3_ 


7-n 


(2.9) 


n 


This expression gives the transfer f 'notion between the 
k-th output and the 1-th inputs For different values of 
k and 1 this describes the traisfcr function matrix Hj^i(Z) 
with each transfer function (Z) as its element. 

To express eqn. (2,9) in the transfer function 
matrix form, wo write it separately for each input variable 
i.e. for each value of index ‘ 1 ', for a given output = The 
total output variable is obtained by using superposition 
theorem . 


P 


P 


"k 


= H “a = H-bi q (2) 
1=1 1=1 


( 2 . 10 ) 


Whiting eqn. (2.10) for different values of ’k‘ transfer 
function matrix (Z) is obtained. 
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2«1«3 A particular cas e 

Ai exanple of a nraltivariable sys’cem having two 
input output variablGs follows ; 


Here 1=1,2 




So, c^^(Z) = 


'k1,o kl,l 


,111 


1 — b Z — b-i^-i p 2 -t**- b. 


— r^(Z);k=1 ,2 


,-1 


and Cj^(Z) 


\2,1 ^ ■*■••••■^ 2.-2 


k 1 

-m 




1-b z'”’ - b Z"^ b. 


-n 


r2(Z) ;k=1 ,2 


'k2,2 ~ k2,n 

These are the output variables due to inputs r^ (Z) and 
r^CZ) respectively. The total output variable is, 

(2) = Gy (Z) + Ci^(Z) ; k = 1,2 

Thus the system is represented in the following transfer 
functic.i matrix form - 


c,^^(Z) +0^2(2) 


C 2 i(Z)-KJ 22 (Z) 


C^(Z) 


Co (2) 


1 




7 -m 


-1 


1 -b. 


H-' 


11 ■m"' ^2.o'*^12'1^ 


-n 


^""^1 2 *1 2”*^ - • • » -b-i 


3^ 

z--n 


^ Ij'^Z" — ••-b^1 “12,1^ "•••"^12,n 

^21 , o’^^21 j >! • +b2^ ^ 0.2 ^22 , o"^ 2,1^ ^ +. . . +a22 

1 — bo-i ^ Z » — bo ^ 2 


-m 


or 


21,1^ 

C (Z) = HCZ) R (Z) 


1 -b 22 j ^ 2 "'^ - . . . -b22 


-n 


'^2 

1 

i: 

:ri(Z) 


1-2(2) 
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2*1 -Jj.en.er.al transfer f-unction m atrix 

As in the secticn (2.1 « 3), the eqn* (2«9) 
may bo used to -write the transfer function matrix of the 
model ha-ving *p' inputs and ’s' outputs. So, the reduced 
order, linear discrete system may he expressed i'n the 
pulse transfer function matrix form s 

C (Z) = H(Z) R(Z), (2.11 ) 

■whero, 

C^(Z) =[c^ C 2 ... (2.12) 

R^(Z) ='^r^ r2 ... rj, (2.13) 


H(Z)=^ 


Wi ** 

-4 ^‘L.u_k,, 7-h \ -L t 


fVv, 

‘-W 


(2o1^^ 


of from eqi. (2.4-), in the form of the equivalent 
difference equation ; 
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Pm P 


k=1,2,ooeS (2.15) 


n p 


j =0 1=1 


\ 

a. T . r + / L b, -, . c, T . 


k = 1 ,2, . . . , s (2o1 6) 


Letting i range from zero to some integer * q' , 
eqn. ( 2 . 16 ) may be represented by the matrix equation : 


V = %,q 


( 2 . 17 ) 


where : 


V \2,0 ••* \p,o’ ^^1,1 


k1 ym 


k « « a, ) (b, ^ i**b, (b ^b, *b, ... 

k2,m kp,m^ k1,1 kP,V k1 ,2 k2 2 kp,2^ 


(b , b, ^ .... b, )1 

k1 ,n 1-52, n’ kp,n^ i 


( 2 . 18 ) 


% „ - Pki - 


(2.19) 
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( 2 . 20 ) 
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Vector ^ is termed ’Parameter’ vector and 
matrix is termed ’Data’ matrix. 

The problem^ then, is the determination of 
(m +n + 1 )p components of the Parameter vector from 
the matrix oqn. ( 2 . 17 ) the left hand side of which contains 
the Data matrix, obtained from the input- output data at 
sampling instants* Its right hand side is the vector of 
total output variables, all of which are known. Therefore, 
the problem simply appears to be that of the solution of 
a set of linear algebraic equations* However, in practice, 
the number of rows in the Data matrix is much larger than 
the number of columns. This is particularly so when the 
data is taken over a large interval and the sampling rate 
is high. Moreover since we want the order of the reduced 
model to be as small as possible, the number of columns 
in the Data raatrix will naturally be small. Because of 
this reason we get more number of equations than are the 
number of unknowns, and the set of linear equations repre- 
sented by ( 2 . 17 ) does not possess a unique solution, unless 
m and n are the actual orders of the polynomials of the 


pulse transfer functions. As this is not the case, by 
hypothesis, we can determine that parameter vector, 

3^ , which gives the ’best’ fit. One suitable criterion 
for an optimum fit is the minimization of the sum of the 
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squares of the errors between the responses of the actual 
system and the reduced model at the sampling instants* 

This is the criterion considered in the present work. 

2.2 Method of solu t ion 

As stated earlier, the problem does not have a 
unique solution. In other words, since the Data matrix 
is rectangular, it does not have an inverse in the ordinary 
sense. However, for 3 > (m+n +1) p, the least-squares 
solution to eqn. (2.1?) ,is given by, 

A LH 

X, = A C (2.21) 

-k \q -q 

LM 

where A^^^ is called ‘Left Minimum Pseudo Inverse ‘ of the 
matrix A^^^ and is defined by, 

LM - LM T -1 „ 

V] ■ 

The properties of the matrix pseudo-inverse have 
been discussed thoroughly by Penrose and GrevilleQ+8l . 

A description of matrix pseudo- inverse is also given in 

A. 

the appendix. There it has been shown that from eqn. 
(2.21) is determined in such a manner that the sum of the 
squares of error between the right hand side of eqn .(2.1 TA 
(which is the measured response of the actual system) and 
its left hand side (which is response of the low order 
model) is minimized. In other words, pulse transfer 
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functions in the transfer function matrix (2.14), with 
the coefficients determinec’’ from (2.1 S), ^reld respon- 
ses of the low order model which approximate those of the 
original system. 

2.3 iU-goritlms based on the matrix pseudo- Inverse 

In practice, the Data matrix is quite large and 

LM 

it is not feasible to calculate directly. This is 
particularly the case when a large amount of data must be 
stored in order to miniraize the mean-square error over 
a sufficiently large interval, especially when the samp- 
ling rate is high. This difficulty can be overcome by 
using the following recursive algorithm, originally 
proposed by Albert and v3ittler l 591 and modified by 
Wells C^O"] for the case where a row is added to 
every time, as happens with each pair of additional input 
output data; 




3 ^ 


and let 


%c^+^ 


(2.25) 


Is, 1+1 


The recursive algorithm that follows is divided in two 
phases* The first phase is the starting phase* 
start from q = 0 and go upto q^(m + n +1) p and 
use the following recursive relations ; 


-k, q +1 ^ 




a Q, a, . 

~k,q+1 h,q“k,q+1 


(2.26) 


where 


:,q %,q+|C\,q %,q-+l] 


(2.27) 


•%s,q+l %,q %,q+1 


^k,q+1 = Pic,q 


Pk,q -^Cjq+i^l^^lSjq ^.q+ll "’"fes ,q%:, q+l] 


-k,q+l '^k,q -^jq+l ^ 

I- k,q 


^»q+il 


'k, q — k q +• 


-i]l\q %s,q+ll C^'^,q+l ^k,q •%s,q+ ll 


. T 


(2.28) 


A 

x^rith the initial conditions X = 0^ = I and P = 0. 

— O “ O o 

The second phase follows the first and starts for q>(m-in+1)p. 
The reciiTsive relations in this phase are ; . 

m ^ 


Y = Y + :_ k.q--^.a+1 k .q^1~ -g 

k q -H -j , p 

^ ^-^c.q+l ^k,q %s,q+1 


XA9-.lir_.-%^_q +1 --kq) 


(2.29) 



3^ 


whore 

P = P 

k,c[-H ^k,q 


[^k,q %,q+l^k,q %,q+l] 
q-M ^k,q %,q+l 


( 2 . 30 ) 


2 «U- Do t o rmin at 1. on of . the e qi3 iwal e n t c ent in -uo-u s- ti me m odel. 

The algorithm given previously determine the 
low-order pulse transfer functioi matrix H(z) for the 
system# The problem now is to determine the corresponding 
continuous-time system transfer function matrix H(t). It 
should be noted that II(s) is not the Laplace transform 
of the impulse response which has the z-traisform H(z). 

The relationship between H(s) and H(z) depends upen 
the form of the input between the sampling instants# 

Anders assumed that the input is held censtant 
between tho sampling instants. This is equivalent to 
the introduction of a zero- order hold, so that H(z) may 
be regarded as the z-transfer functioi of H(s) preceded 
by a sampler and a zero-order hold# If the input to 
the system is a unit step this does not introduce any 
errors at all. 

Smith £613 has given a method for estimating 
the Laplace- transform coefficients of a dynamic system 
from sampled-data values of its input and output. 

Aivantage of this method is that it is quite straight 
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forward and can be easily programmed in a computer as 
tho algorithm is iterative in natureo The method is 
applicable to tho systems with roise. Its application, 
on ndso-free systems is describ.d in section 2«4.1. 


Description of Smiths s m etha~l 

In this method the dynamic system is assumed 
to he linear vilth a Laplace transform of the form, 


H(s) = 




+ 


^-] s + 


+ oCp s 


P 


^ 0 + s + ..o + [3 s°- 


(2.5+7) 


The coefficients expand 

in this transfer function are not h.novm. The pulse 
transfer function H(z) which exactly relates the input 
and output sampled data sequences is 


/s 

H(z) 


A 


an + a- z’ 




-m 


/V A **1 A 

b + b. z + ... + b z 
0 1 “ 


(2.4-8) 


In general each of tho cooffscient in this 


pulse transfer function is a function of coefficients 
of the transfer function (3.4-7) . In other xfords. 
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/s /\ 

^•O " % ^ ^ o» ^ 1 »• ••> p Po’ (3q) 


A /s. 


m m 

A A 


\ p’ Po’ Pq) 


^o ^0 p’ Po» Pl»***’ Pq^ 


(2.if9) 


A A 

^n ~ \ ^^Q»'^-|?***’'=<p’ Po’ Pi’***’ Pq^ 

Noxif the pulse transfer function (2.4^) must be 
of the same form as the pulse transfer function HCz), 
estimated from the data# If it is not, then the form 

of H(s) is changed in such a manner that the forms of 

/\ 

the nulsG transfer functions H(z) and Hfz) are exactly 
similar. 


The pulse transfer functicn H(z) ostimated 


from data is given by, 


H(z) 


Oq "t" q*! z ^» ♦ • + a^Q z 
bq "^^1 2 "'^+ z ~^ 


(2.50) 


The coefficients, a^, a^,..o, a^ and b^, b^,..., 
b^ of this transfer functicn are knowpo 

Now equating the coefficients of H(z) and H(z) 
for corresponding powers of Z, we get a system of 
equatiens, withcx^jc/i ,.«*Mp Po’ Pi?*’*’ Pq 
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unkn owi s , 


% 

<b< o’^1 ’* *• ’^P’ P 

0» ?••• » p 

II 

A 

o’ ^ p’ 



b 

0 

• 

Po’ Pl ’ • •• ? 


« 

• 

A 




^n 

(o<Q, o( -]>•••> o<p? 

po’ 

Pq) = h 


Solving eQns* (2. 51) for the untoowns, oC 

o I 

and p^,..., P g ^5 the coefficients of the 

transfoi' function H(s) are obtained. An example of 
this method is given in section 2. 4-. 2. 

2.4.2 1 cxajiiple of the apt lication of Sml ~ ^h* s metho d 

Suppose that the pulse transfer functicn 
obtained from sanpled data is, 
b^ +b^ +lD^ 2“^ 

H(z) = (2.52) 

1 + a^ 2 + ap 2“ 

Here, the values of coefficients b^, b^ , b^, a^ and 
32 dre known* The form of this transfer function 
suggests that H(s) to bo estimated from it should have 
the form m , 



39 


oLo s 



( 2 . 53 ) 


( 2 . 5 ^) 



( 2 . 55 ) 

( 2 . 56 ) 

( 2 . 57 ) 



As can be seen, there are four unknowns, 
and ^2 equations relating these 

variables are available. Thus exact values of all of 
these coefficients cannot be determined. To get approximate 



^0 


values, 
Step I 
Step II 

Step III 

Step IV 

Step V . 


of following iterative scheme may be used, 

- Findp^ from eqn. (2.58) 

- Assume a value of and find-t^^ from 
eq.n. (2.;,^), 

- Combine eqns. (2.56) and 2.57) to eliminate 

trign ometrical expressions and by substituting 
the values of , (5 ^ and find'X^. 

- Substitute values of< 5 <^, in eqn.(2.55) 

and get a value ofoc • 

o 

- Compare the values oft^^ obtained in Step II 
and step IV, If they are close enough, stop, 
otherwise, assume a new value of f3^ and 
repeat the above steps. 

A computer program can be made to inplement 


the above scheme 



CHAPTER III 


■fiP PL I CATION TO A NUCLEAR REACT 0^1 TURBINE MODEL 

3,1 B rief d e.soriTDti on nf the_nucl.par re act or turbine model 
Tlie development of the model of the nuclear 
reactor turbine given here is based on IBM report £62“\ and 
reference 1631 . 

The block diagram representation of the turbine 

$ 

coupled to the nuclear reactor with its various elements 
is shoim. in Tig. 3.1 • The steam turbine assembly can be 
divided into various subsystems like high pressure, 
intermediate pressure, low pressure units, feed water 
heaters, etc. The well known process representation for 
this part of the modified Rankine cycle is represented in 
Fig. 3*2. The steam flowing from the reactor through 
the main steam valve is just dry, hence the steam exiting 
from the high pressure turbine will be wet, as indicated 
by the condition of steam at point 2 of Fig, 3*2. The 
wa+-or part of tlw) wet steam is removed in the moisture 
separator and the dry steam at pressure pj^ is superheated 
to temperature % in the reheater by using secondary steam 
drawn from the main supply line through a Secondary Steam 
Valve. This superheated steam expands ttu’ough interme- 
diate pressure and low pressure turbines and then enter'’ 
the condenser. The condensed steam, ice-, feed water, 



flov^s through bhe feed water heaters before it enters the 
steam Crum* The feed pump ■’’n between the tv/o heaters 
raises the pressure from the condenser pressure to the 
reactor drum pressure. For the above system the model 
and plant dynamics are described below. 

3*2 Model an d the Plant dynamics 

The model was developed by the application of 
well known physical laws and a few pertinent empirical 
relationships* Various assumptions made in the develop- 
ment of bhe model are presented first, followed by the 
model itself. 

(1 ) Tlie relationship between enthalpy, specific 
volume and pressure (i.e. the equation of state) for the 
superheated steam follows Callender's empirical relaticn- 
ship for all noi-ccnstant pressure processes. 

(2) The I.P. and L.P. turbine dynamics are 
lumped together, as both have a correspondingly slower 
steaiT! flow dynamics as compared to the H.P' turbine. 

(3) The flow through the throttle valve is 
assumed to depend on the square root of the upstream 
pressure and density of the steam and linearly on the 
area of flow. 

(k-) The regenerative bleed flow depends linearly 
on flow rate. The flow-time relationship for the turbines 
is represented by a first order lag. 
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( 5 ) In order to take into account the reduction 

of rotabional losses 5 root and tip clearance losses for 

the TI.P. and L.P« turbines at reduced flow rates, 

V} and H 5 are nialtiplied by correction fact cars 
•flp Lp 

and i- , respectively, which are assumed to 
^ CFLP 

depend linearly on the flow rate* 

( 6 ) The bleed flow is tapped right at the end 
of the high prossuro turbine, and thus the entire flow 
through the turbiiG participates in p reducing torque 
while only one-half of the bleed flow is assumed to 
produce torque in the L.P. turbine. The secend assumption 
is duo to bhc fact bbat bleed f].ow is removed at various 
points in the I.P. and L.P. turbines. 

(7) Heat exchange in the reheater is assumed to 
bo perfect and the dynamics of mass balance and energy 
balance are lumped at a single point. Pressure drop in 
the reheat er is assumed to be zero. 

Hdth these assumptions the system equatiois 
representing the turbine dynamics are given below. The 
variable names are explained in Nomenclature B. 

(i) Nozzle Chest 


w-,hg - W2h^ ^ ^ 


^ V 
c c 


772 dt' 


1- k J 
1 


(3.1) 



~ ^2 

dt V 

c 


w 


1 


c* A p 
1 1 




A* (p^ Pn - PaPa) 


X = (hg - Pa)/bfg 

P = 1 

^2 x" (1-x) Vf 

V|. = 0 * 01 84 

Vg = 2,288 - 0.0166 (p^ - 200) 

= 12.57 X ro^ - 27 (p - 1000) 
Po =/o ^6 “"o - V 


(ii) I-iLgli pressure turbine 

^^2 ar " ■ "bh ^ - ''2 


Wv 


Blip 


BtP ^2 


(iii) Moisture separator 


^MS 


! 



<'^2 - ShP '' 2 ^ - "2 

^2 " " 

— 

hfg 2 

355 + 0.4+ (Pp^ - 200) 


X 1 .055 X 1 0 
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(3.2) 

(3.3) 
(3 *4) 

(3.5) 

(3.6) 

(3.7) 

(3.8) 

(3.9) 

( 3 . 10 ) 

(3.11) 


(3^1 2) 
(3o13) 


(3.14) 
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- 0.4- (p^ - 200) X 1 .0^5 X 10 


5 


(3.15) 


(iv) Re heater 

(a) Main steam 


R 




dt 

dt 


Y-r 


= [^ 


+ hg - W3 


P r 

■^R R 

^R “ ^R \ " ^2^ 
1 /2 


Wo = k. 


R^ 


(w^ - w,) 




R R 


(3 .16) 


1-K^J 


( 3 . 18 ) 

( 3 . 19 ) 


h = [1198.4- + O.Of (pj^ - 200)^ X 1 .055 X 10^ (3.20) 


(b) Reheater steam 


^PR _ WpR-wlj, 

- ^7- 


(3.21 ) 

<1 Or _ ''pR + 7 , , 

dt - 2 T iT - %) 

^ ■^R2 % 

k. 

^R2 

(3.22) 

Wp^= p 


( 3 . 23 ) 

Tg = 544-.6 + 0.125 (p - 1000) 

+ 1+60 

(3.24-) 

^R ^ ^ 


(3 .25) 

R = 201,7/0.365 (520 +4-60) 


(3.26) 
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(v) I 

•P • and L. P 

. turbine 


dw^ 

dt 

: (1 „ p ) 

1 

(3.27) 

""blp 

“ ^BLP ''^3 


(3.28) 

(vi) 

Heater 1 



diE 

_ ^1 

- Hi 

p 0 fw 

T 

•^H1 

(3.29) 

Qhi = 

He (w* + 

^ lip 

''blp’ 

(3. 30) 

(vii) 

Heater 2 




_ ^Ez 


(3.31) 

dt 

%2 ^Fw 

^H2 

If 



(3.32) 

dw^ 

ar = 

t 

^ "^PR " ^HP 

(3.33) 

^H2p 

^H2p 


(Ttii) kvork output equations 


'bp = '^2 (‘‘o - h’ ) 

% 

(3.34) 

LP ~ ^3 

(3 .35) 


=]j067 + 0.37 (% - 200) - O.OoH (P:^- 200)^- . 

O 

0.1 (p^ -1000^x1.055x10 (3.36) 

h}^ = 898 X 1 .055 X 1o3 
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* 


HP -^GFHP HP 

(3*37) 

H LP ~ GFLP ^ LP 

(3.38) 

n*=n =0.86 

LP HP 

_ 

(Wj/Wg) ) 

wUw* - 

4.CFLP = / , ' ^ 

(v^/w^) (1- w^) 

( 3 . 39 ) 

(3.4-0) 

= 778/(1 .075 X 1o3) 

(3.4-1) 

TORQUE = T + T 

HP LP 

(3.^2) 

TORQUE X SI 

POWER = 

500 

(3.^3) 


A closer look at the equations reveal that 
equations (3.1), (3.2), (3.10), (3.16), (3.1?), (3.21), 
(3*22), ( 3 . 27 ), ( 3 . 29 ), (3«"1) and (3.33) "xe differential 
equations with the rest being algebraic relationships. 
Thus, in the turbine model there are 11 nonlinear diffe- 
rential equations and 32 algebraic relationships. 

They can be represented in the state vari^le 

form as, 

i = 2) (3,Mf) 

where, 



T 

statevector , X = (X^ , ... X^^) 

= bfj 14^ Wgp 

(3A5) 

input vector, = (U^ U2) 

= (Aj (3.^6) 

3.3 Low order molel of the plant system 

The method developed in the present work w as 
applied to obtain lav/ order model for the nuclear reactor 
turbine system described above. ALthongh^. for this system 
the model equations were available, in actual practice, 
there is no necessity for these equations to be present 
in order to obtain low order model for the system. Measure- 
ment of system responses at discrete intervals of time due 
to known excitations to the inputs is sufficient to deter- 
mine low order model for th^ system. In the present work;, 
this data v/as generated from the model equations. 

As can be seen from the equations describing the 
plant dynamics of the system, it is not easy to find out 
the exact order of a transfer function relating an output 
variable to a particular input variable. However, looking 
at the number of differential equations in the svstem, 
which is eleven, it can be assumed that no element in the 
transfer function matrix can be of order more than eleven. 
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It was decided to model each element in the transfdc 
fijncticn matrix as of second order and go for higher 
orders if the respoises of the model did not match with 
the responses of the original system. However, strategy 
was to go for even lower order model and find out the 
lowest order model possible. 

In practical situations if plant system equa- 
tions are not available, one can decide on the order of 
the transfer functions and aiopt the above procedure to 
obtain low order model. 

Having decided to model each element in the 


transfer function matrix (2.14) as of second order, the 
most general second order transfer function with m = 2, 
h = 2 in eqn. (2*9) was selected. Thus a typical 


second order transfer function will look like, 

H. ^kl,1 2"’^ -^^1,2 2“^ 

H ^ i ( 2 ) - -—. 1 - -*2 

^ ~ \l,. 2 - V,2 2 


A 3*3 percent step change in the main throttle 
valve area (A^) and a 10 per cent step change in by-pass 
valve area vere made for response studies. Changes 
over the steady s tat e valves in the output variables dne 
to both of these inputs acting together, as well as 
individually, were taken with the sampling interval 
T = 0*1 sec. and using 200 samples* 
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These output values were put in the Data matrix 
(2.20) and equatican (2.21) was solved f or unl?nown 
coefficienbs hy using the algorithms given ±a section 
2.3. The computations were carried out on DEC 1090 
Computer System. 

The coerficients of the pulse transfer functions 
were converted into corresponding centinuous transfer 
function coefficients hy the method described in section 
2 .tf- • 1 • 

Other low order models for the s'^stem with 
transfer function matrix having all first order transfer 
functions as its elements, and with transfer function 
matrix carrying a combination of first and second order 
transfer functions were also derived. Discussion of 
responses of these various models is givai in Chapter IV. 
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DISGUSSIOH ON THE RBSCTLTS MD SELECTION OF 
LOW CRDER. MODEL 


^«1 Intro d uctio n 

The present method was applied to obtain low 
order models of the nuclear reactor turbine system 
described in Chaptoc III* ]h t^^e present chapter 
different low order models for the system are presented. 
R.esponses have been obtained for - 

(1 ) Reheater temperature due to 3»3% step change in 
the main throttle valve area and ^ofo step change in by- 
pass valve area. 

(2) Low pressure turbine torque due to 3*3% step change 
in the main throttle valve area and 10^ step change in 
by-pass valve area. 

(3) High pressure turbine torque due to 3.3^ step diange 
in the main throttle valve area and ^0% step change in 
by-pass valve area. 

(^) Steam flow rate due to 3*3^ step change in the main 
throttle valve area. 

3*3^ and 10^ step inputs are chosen only for the 
the sake of comparison of responses cf the low order model 
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and the acbual system, the responses of which are known 
for these inputs 1623. The methoi at such does not 
put any restriction cn the input. 

4.1.1 Heheater temperatur e 

The four models obtained for this response 

are - 

Model I - Both main and by-pass throttle valveshaving 
second order transfer functions. 

Model II - Main throttle and by-pass valveshaving II and 

I order transfer functions respectively. 

Model III- Main throttle and by-pass valves having I and 

II order transfer functions respectively. 
Model ly - Both main throttle by-pass valveg having I 

order transfer functions. 

These models are presented in table 4.1 (a) 
and their responses are shown in Fig. 4.1 (a). In this 
figure, it can be seen that the best over^l response is 
obtained by model I. However, the best peak response 
is obtained by model II. This is because of the ncn- 
oscillatory nature of the first order transfer function 
(used to mod^ the behaviour of by-pass valve) super- 
imposition of which tends to reduce the overall peak 
overshoot. The transient responses of Models III and. IV 
are considerably in error, lykgnitude of steady state 
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error (.07^) is same in all the models. 

Tabic ^+.3 gives time constciits (I order trens. 
f-unc.) •’nd settling timo (II order tr^ns* func.). Talsle 
^*4 givjs the CPU time for model generation. 

Computer Processing Unit (CPU) time for the 
various models are shown in Table 4.5+. 

t|-«1 «2 Low pressure turbine torque 

The models obtained for this response are : 
l-lodel I - Both main throttle and by-pass valveshaving 
sec Old order transfer functions. 

Model II - Transfer functiens of main throttle and 
by-pass valves of second and first order 
respectively. 

Model III- Transfer functions of main throttle and by-pass 
valves of first and second order respectively. 
Model IV - Both main throttle and by-pass valveshaving I 
order transfer functions. 

These models are presented in Table ^.1 (b) 
and their responses are shewn in Fig. 4.1(b). Ag can be 
seen from this figurcj the best transient and steady- 
state (.07^ error) responses are obtained hy Model I. 
Response of Model II gives much higher steady-state 
error (.21^). The transient response is also not as 
good. Much poorer transient responses are obtained by 
Models III and IV. However, the latter gives much 
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smaller steady-state error (.25?? in model IV compared 

to •42^ in Model III). Time constants, settling time and GPU time 
c’.re '■'■-'.ven in Table 4.3 and 4.4 respectively. 
lilgd pressure turbine torque 

Foitr models are obtained for this response and 
they are ; 

Model I - Both main throttle and by-pass valveshaving 
second order transfer functions. 

Nbdcl II - Transfer functions of main throttle and 
by-pass valves of second and first order 
respectively. 

Model III- Transfer functiens of main throttle and 
by-pass valves of first and second order 
respectively. 

Model IV - Both main throttle and by-pass valveshaving I 
order transfer functions. 

These models are given in Table 4.1 (c) and 
their responses are shown in fig. 4.1(c). The figure 
shows that the best steady- state as well as transient 
responses are obtained by >fodel I. But the best peak 
response is given by Model II. This is because the 
peak overshoot is reduced by the ncn-oscillatcry behaviour 
of the first order transfer function which is used to 
model by-pass valve characteristics. Model III gives 
higher transient and steady- state (.65??) errors, than 
the other two models. However, the highest steady-state 
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and fcrsnsient response errcrs ere given by Model IV. 

Time constants, settling time and CPU time are given 
in Table ^*3 rospectively* 

Steam flow rate through nozzle chest 

Two models are considered for tbis response ; 
Model I - Second order transfer function 
Model II - First order transfer function. 

Those models are given in Table 4.2 and their 
responses are shown in Fig. 4.2. As can be seen from 
the figure, Model I gives the best transient and steady- 
state responses. Model II gives much more transient ecror, 
however, steady-state responses of these two models 
do not differ much (.10^ error in Model I and .15^ 
error in ^bdel II). 

TeOole 4.3 gives time constants' (I ojder -trang* 
fLii.c.) 'and settling tidg time (II order trars. func.). T 
4.4- gives the CPU tine _ for models generation. 

4^,2 Se .le ction of . a l ow. , order model 

In the previous section several low order 
models fer the nuclear reactor t'-'rbine system were 
obtained. In ti is section suggestions have been male 
regarding the selection of a particular model. 

The low order model is selected on the basis 
of application which it is expected to be put to. In 
the on-line applications like adaptive control of a 
system, the need is to generate in smallest possible 



time a '' ow order mcd^L, which will d.^-scribe the responses 
of the original system with sufficient accuracy* Here, 
the tjjuo i'Oqaired to generate the low- order model is 
the most inportaii: factor. Thus a model with smaller 
number of parameters will be given preference over a 
model \'/ith more parameters even irhsugh its response may 
not be as accurate as that of the more complex model. 
However, in the situ.ations where the need is just to 
obtain a low order model on the system, and to continuously 
refine it as more information on the system becomes 
avaiJ.able in the form of input-output data, the time 
required to generate the model is not a very important 
factor. In these cases, the low order model which 
gives the most accurate responses will be selected. 

Based on these considerations we will now select 
low ordor models of the nuclear reactor turbine system. 

V .2 .1 Most anpr onriate mode l for on-line app licatl m^. 

As can be seen from Table generation time 

of Model II for all the four rcspon'ses - reheater tenpe- 
raturo, low pressure turbine torque, high pressure 
turbine torque and steam flow rate through nozzle chest, 
is around 60^ of that required to generate Model I. 

Thus for all the responses, Model II is selected and 
the low order model of the system for these responses 
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is, 


1 


‘ 9.9896 (1 + 0.9681 s) 

8<,2112 


s^ J 0 . 6624 -S + O.asc’"- 

1506 .1 1 d: (1 + 0.6605 

s + 0O83V ■ 

s) 99-8512 


s^-^. 662 i+ s-K}, 23 &i- 

s + 0 0 0^ 1 8 

1 

■c ] 


2511 • 87 d+ (1 +0.6831 s) 

1 0.03+ 

3 1 

^ ] 


■ s^ -HD. 6562 s+ 1 .2061 

5*3287 

s+O.08+0 

sHb*0830 ■ 

0 


^2 




where, 




Reheater temperature over its steady 
state value 

Low pressure turbine tcrque over its 
steady state value 

I'E-gh pressure turbine tcrqu'^ over its 
steady state value 

Steam flow rate through tho uozzle chest 
over its steady state value 
Main throttle valve area 
By-pass valve area* 
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4«2»2 >fost appropriate model far off-lin e applicatioi. a 
Fcr off-line appli aticns, time r''iuire(i to 
gene3-'abe a model is not the most important factor, llhe 
model of reasonable order which best describes the 
system respoises are selected. So, for the nuclear 
reactor turbine system. Model I is selected for all 
the responses. Hence the most suitable model for 
this system is, 


1 


’'9»0225(1 + 0.861^s) 

8.7032(1 -HD. 1297s) 



s^+ 0.69!;2s +0.29^9 

?+ 0.6830s+Oo2315 , 


^2 


1 521 .7381 (1 +0,^32s) 
s^-K). 6222s +0.26 01 

1 00.6261 (1 +0o5-1 18s) ; 

«.625-8s+OJ936 i 

,>i'‘ 



21^99.085-2(1 +0. 5369s) 

11 .5255- (+ +0,671 3s) 


S' 


+0.6721 s+ 1.2555- 

s^-i0.63^s+O,8e)+9 

i k 2 

. . 


5.5-1 ^(1 +0,29673) 
s^+0. 6971s +0.7225 

0 

i 




IV 


8. 979*+ 
3+0.0819 


8.1 312 
s +0.0806 
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Table 4.1 (b) 

Low Pr e ssure Turbine Torque 


Input 5 

^'^'vs.ijra^able | 
Model \ 

Change in the inaln 1 Change in the by-pass 

throttle valve area 1 valve area 

i 

t 

I 

1 521 .7381 (1 40.4432s) 

40.6222s 40.2601 

100.6261 (1 40.4118 s) 
s^ 40 . 6248 s 40 . 1 9 36 

II 

1 506*11 C4( 140.6605 s) 

99.8512 

2 

s’ 40.6624 s 40. 23 c4 

s 40 0 08 I 8 

III 

2009.8716 

s 40.0809 

98.731 2(1 40.31 21 s) 

3^40, 621 2s 40. 2025 

IV 

L 

200 7. 081 3 

S 40 . 0826 

99.2213 

: s+ 0,0832 
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Table 4.1 (c) 



Turbine Torq ue 




















Table 4.2 

Steam Flow Rate Thrcfugb No?;zle Chest 


64 



II 


5.3287 

s + 0. OS+0 
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Table 4.3 

Settling time (II order trans. func.) end Time 
Constants (I order trans. 


ilespon-O ] 

\se 5 R.T. 1 

Model^ — — t J 

1 

L.P.T.T. 1 

H.P.T.T.^ 

S.F.R. 


II 






order 

1 1 ikhr 8 

12.85 

11.98 

11.50 

I 






JL 

II 


' ' — 1 




order 

1 1.71 

1 2.80 

12.65 



II 


i — " — ^ — -n 



II 

order 

12.08 

12.1 7 

1 2.10 

11.90 


I 




■ - 


order 

1 1 .99 

1 2.22 

12.20 



m 


. H, . .. , , 


- 

lU 


12.23 

1 2.36 

12.18 

- 


II 



■■ 

- 


order 

1 2.40 

12.88 ; 

B 



I 





IV 

order 

12.20 

1 2 .1 0 

12.1 9 



I 




- 


order 

1 2.51 

1 2.02 

12.30 



j&bre-viaticns ; i4T. - Reheater Tender ature 

L.P.T.T.- Low Pressure Turbine Torque 
H.P.T.T.- High Pressure Turbine Torque 
S.F.R* - Steam Flow Rate# 
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Compirber Proemsing Unit !B.me (sees.) 


Input 

^\variablo 

r f ■ ' ‘u '1 

H.I. ( L.P.T.T.S H.P.T.T. 1 S.P.E. 

5 5 5 

1 i 5 __ 

I 

92.13 

91 .18 

92.22 

60*38 

11 

55.25 

55.51 

55.61 

36 .^-3 

III 

55.90 

i 

55.72 

I 

55.11 

- 

IV 

i 

35.28 ! 

I 

35.09 

35.11 



Abbr GPiaticn s : R.T. 


L.P.T.T. 

II.P.T.T. 


Rehoator Temperature 
Lov Pressure Turbine Torque 
Pligh Pressure Turbine Torque 
Steam Plow Rate. 


S.F.R. 
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C ONCLUSION^ LIMITATIONS Off THE PRESMT METHOD MD 
SCOPE FOR FUTURE WCBK 


5»1 Goncluslcn 

In the present -wcsrk, a method has “been developed 
to obtain low order models of large multi variable systems. 
The method is applicable to both linear as well as non- 
linear s/stems end does not require the transfer function 
or vector differential equation of the system to be 
specified. Only measured input-output data of the system 
at discrete intervals of time is required to determine 
its linear, low order, discrete- time model. The input- 
output data is put in the form of a svstem of difference 
equations, with model parameters as uh’^nov-jns. The equa- 
tions are solved for these unknowns by the use of matrix 
pseudo-inverse • The parameters thus determined, minimize 
the sum of the squares of the errors between the responses 
of the original system and its low order model at the 
sampling instants. A recursive algorithm is used for 
the solution of the system of difference equations hy 
matrix pseudo-inverse technique. This reduces the 
computational and data storage requirements and enables 
successively refined model of the system to be determined 
as new data from the system becomes available. 
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in iterative technique developed hy Smith |I 59 l 
is usei to obtain the corresponding continuous- time 
model of the discrete-time model obtained by the above 
method. 

The method was applied to a nuclear reactor 
turbine system £623, ^633* which has 11 non-linear 
differential equations. Reheater temperature, low 
pressure turbine torque and higlh pressure turbine torque 
responses for 3 * 3 %* step change in main throttle valve 
area and 10 ^ step charge in by-pass valve area and steam 
flow rate through nozzle for 3 * 3 ^ step charge in main 
throttle valve area were used to obtain the low order 
models for the system for these responses. The different 
models, with both transfer functions of II order, both 
transfer fnncticn of I order as 'vrell as conbinaticn of II 
and I order transfer functiois were obtained. Responses 
of these models have been discussed and particular models 
most suitable for on-line and off-line applications have 
been suggested. 

5*2 Limitations of the qresent method 

The limitations of this method are as follows, 

( 1 ) Steady state error between the model and system 
responses cannot be predicted before obtaining the model. 
This is because the method minimizes the sum of squares 
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cf the errors between the system and model responses 
and ■" fixed expression for steady state error, like 
the one given in Davison's method^ 2^3 cannot be 
derived • 

( 2 ) Least squares criterion may not giveil the best 
results in all systems* 

( 3 ) Values obtained from the iterative algorithm may 
not converge easily. 

ih) Bie time constants of all the transfer functions 
may not be the same. Thus order of the low order model 
may not be what was specified. 

5»3 Scope for future work 

Several suggestions are given below which may 
make the present method more useful j 

( 1 ) In the present work, the method was applied to 

a nuclear reactor turbine system. It should be used to 
obtain low order models of several other systems to give 
more confidence on its applicability. 

(2) The method may be extended to be able to deal with 
systems -whose parameters vary slowly with time. 



(3) A t package’ computer program can "be developed xdiich 
should directly give low order ^ continuous time model 
hy doing the lov;' order modelling and converting the 
discrete- time model into the corresponding continuous- 
time model. This program shoild also select the most 
appropriate low order model for a given application 
by considering the response errors and model generation 
time. 

(A) The methoi may be modified for 'control law* 
reduction applications. 
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APPENDIX 


!IA!IBIX P lEaPO-INVERSE 


1 » Concept of llatrix Pserao Irfverse 

For a rectangular vjatrix, an inverse does not 
exist in the ordinaiy sense. It is impossible to have 
a unique and exact solution for a set of siimiltaneous 
equaticns when the number of equations are not equal to 
the number of unto owns. Tv/o cases might arise ; 


Case (a) t Number of equations exceeds the number of 
unto owns* Consider the equation. 


H X = y (n > m) 

*fMr 

(nxm) (mx1 ) (nx1 ) 

where, 

H is a (nxm) matrix 
is a (ffix1 ) vector 
and y is a (nx1 ) vector. 




H = 

- 1” 

j 7 = 

■ 10 ' 


: 2 . 


4- 


( A-2) 


then eqn.(i'—2) implies either 
X = 1 0, or X = 2 , 
and we cannot find a solution for x* 



The problem, then, is to find a (mx1 ) vector 
X, which approximately satisfies eqn.(i'^1 )• The solu- 
tion is given by. 


x» = X (^3) 

L 

H is a non-vnique (mxn) matrix called Left 
Inverse of H such that. 

H = I21 (A-if) 

where, 

Im is an (mxm) identity matrix* 
letting h^ h^J^ for the system of eqns* (2) we 

find that, 

h^ +ai2 = 1 (a-5) 

Therefore # is given by all points cn the line H-H in 
Fig. A-1 ; hence is not unique. The vector His 
shown lying along the line P. In other words, ai infinite 
number of x's exist, none of which are solutions of 
the original problem, 

X’ = 2; 

However, there exists a unique solution for the minimum 

LM 

Left intense, H given as, 

A ^LM 

S = H y 


jfM 


exists when H has rank m or 


H 


•LM _ 


= (H H) 


-.1 T 

' H 


greater, and is given by, 

(A-7) 
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We will show that the distance between Hx and 
y is minimized for the "best” solution x i 


(jHx - yll > i| H X - yj{ 

As an example, take systems of eqno (A-2), 


(A.8) 


here , 


T 2 2 

H H = (1) + (2) =5 


H™ 

A 

X 


= (0.2) C1 

= Co. 2 




H X for this case is givai by point p^ in Fig» A-1 . 

Pq amcng all points on pp, is closest to the tip of y. 
lies at the intersection of the lines HH and PP, 

To prome inequality (A-8) , we let x be any 

/V 

solution and x be the minimal norm solution, and start 
with the equality, 


- yjj = II H (x - X) + Hx - y)| 

= l|H (x - ^ 2:1! 

= ||H (X - X) + y 11 + 

2Ch(x-x) 3’^ (EH^h ;-^)y (^,9) 

The last term of right-hand side of eqn, (.A-9) vanishes; 
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T ^ ^ 

£h(s - x)3 (HE™ - 1^) Z = (X - X) H (HH^“ - 4)Z 


A T T LM 

= (x ^ x) (H HE - H )2; 

A T r- T' T» 

= (s - X)- Lh- - n^-^z 


= Oo» 

Therefore, eqn. (A- 9) reduces to, 

Ite - 11} =l|H(x - X) N +\\m'z-z\\ 

= l|H(x-x)i( +llHx-Z}{ (A-10) 

Since the first term on the right-hand side of this 
equation is non-negative, the ineqaality (A-8) holds* 

Case (h) : The number of unlaiovms, exceeds the number 
of laiowns* Coisider the system of equaticxis , 

H x = 21 (n>m) (A.-11 ) 

(mxn) (nx1 ) (mxl ) 

We cannot solve for x, because, the determinant 

Ih} = 0. The problem, then, is to find a (nxl ) vector 

X such that eqn • (A-11) is satisfied* 

fhe solution to this problem is given by, 

X = 21 (iU12) 

The (nxm) matrix I? in this expression is called the 

R 

■Eight Inverse of H* H exists ^\rhen the rank of H is m 
or greater and is determined by. 
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1 

■m 


■where, L is the (mxn') identity matrix* 

TD 

As in the case (a), H‘ and hence the soluticn 
(A-.12) is not uniq-ae. There exists, however, a unique 
solution, for v;hich tl:e norm of x is minimal. Such 
a solution x , called the Idnimal Solution, is determined 

- r 

by the Mnlmura Ught Inverse, so that, 

^ IM 

x = (iui3) 


■where , 

H™ = ( H 






